Spin Waves in Disordered III-V Diluted Magnetic Semiconductors 
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We propose a new scheme for numerically computing collective-mode spectra for large-size sys- 
tems, using a reformulation of the Random Phase Approximation. In this study, we apply this 
method to investigate the spectrum and nature of the spin-waves of a (III,Mn)V Diluted Magnetic 
Semiconductor. We use an impurity band picture to describe the interaction of the charge carriers 
with the local Mn spins. The spin-wave spectrum is shown to depend sensitively on the positional 
disorder of the Mn atoms inside the host semiconductor. Both localized and extended spin- wave 
modes are found. Unusual spin and charge transport is implied. 



I. INTRODUCTION 

Diluted Magnetic Semiconductors (DMS) based on III- 
V semiconductors doped with Mn have attracted a lot 
of interest recently, after critical temperatures for the 
onset of ferromagnetism of the order ^i-iUO K have been 
found in Gai-j^Mn^^^As, for x = O.OSS.yBB More recently, 
critical temperatures larger than room temperatures have 
been reported in Mn-doped GaN, enhancing the hope fca; 
extensive technological applications of these materialsJj'El 

While there is general agreement that the magnetism is 
due to charge-carrier mediated, effectively ferromagnetic 
interactions between the Mn spins, there are various the- 
oretical models attempting to understand their detailed 
behavior. The DMS are alloy systems, with inherent 
positional disorder of the Mn atoms; further, spin-orbit 
coupling may play a significant role for hole doping. A 
theoretical treatment which takes into account all these 
factors, and their effects on the magnetic and transport 
properties of DMS, is not yet available. Instead, the- 
oretical models have tended to concentrate on different 
aspects of the problem. 

One class of models, where much work has been done, 
neglects both the disorder and spin-orbit coupHng effects. 
At large carrier densities, where Coulomb potentials of 
the impurity (Mn) atoms are effectively screened out, and 
other disorder effects (such as those coming from As anti- 
site defects which are believed to be the cause of the large 
compensation observed experimentally) can be neglected, 
the former at least can be justified. In such a case, holes 
occupy a Fermi sea at the ton-pf the valence band.la S|tjud- 
ies of the spin-wave spectraQ, Monte-CarlO|-studiesB as 
well as dynamical mean-field theory studicsEI have been 
performed. Overall, their results are rather similar to 
the physics one would expect to hold in a conventional 
ferromagnet. More recently, it has been proposed that 
spin-orbit coupling of the valence band, in the high car- 
rier density (strongly metallic) regime leads to an RKKY 
coupling between Mn moments that is anisotropic in spin 
space, and where the positional disorder effectively leads 
to random anisotropy.liiJ This would lead to frustration 
effects on the magnetic ordering, and therefore unusual 
ferromagnetism. A number of very useful ab-initio stud- 



ies have also been published.0 

In parallel, we have studied effects of positional dis- 
order on the magnetic ordering (in the abse*^p©£-S|pin- 
orbit coupling) of an impurity band modeljl3llj'ESE3 de- 
veloped from previous work doip,™ the context of in- 
sulating II- VI DMS compounds .113113 Such an approach 
should be of relevance at low doping concentrations x, 
at and below the metal-insulator transition, and possi- 
bly even above. Evidence for the existence, of impurity- 
like states is provided by ab-initio studies ,1111 which found 
that occupied hole states near the Fermi level have wave- 
functions mostly concentrated on and near the Mn im- 
purity. More recently. Angle Resolved Photoemission 
Spectroscopy has revealed the existence of the impurity 
band in Gai_:rMna;As withj.^ — 0.035 (very close to the 
metal-insulator[-transition).ll3 A scanning tunneling mi- 
croscope studytj demonstrates the existence of an im- 
purity band in (Ga,Mp)As samples with x=0. 005-0. 06. 
Optical spectroscopycJ also identifies the impurity band 
for x^O.OOOl and x=0.05 samples. 

Electron densities in localized states, as well as states 
close to the mobility edge, can be far from homogeneous, 
unlike Bloch waves. At low Mn concentrations and high 
degree of compensation (which comes from charged cen- 
ters) seen in DMS, short length-scale density fluctuations 
could be significant. This can lead to inhomogeneities in 
the local charge densities at different Mn sites, which in 
turn leads to (microscopically) inhomogeneous magneti- 
zations of the magnetic ions in the ordered phase. Such 
inhomogeneity would be expected to alter the nature of 
the collective magnetic excitations (spin waves) of the 
magnetically ordered system, which in turn would affect 
charge transport through the magnetic excitation pro- 
cesses which give rise to spin-flip scattering. 

In this paper, we develop a numerical scheme to cal- 
culate the spin wave spectrum for a finite but large size 
model of coupled fermions and spins, in the presence of 
quenched disorder. We show that the method accurately 
reproduces the results for a lattice model of small size 
obtained from a standard treatmeatipfi-spin waves via 
the Random Phase Approximation.ciESS We then ap- 
ply the method to the model of Ref. and study in 
detail its collective magnetic excitations. 

The plan of the paper is as follows. The model Hamil- 
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tonian for DMS is described in section II. Section III is 
devoted to developing the numerical scheme. The accu- 
racy of the scheme is tested by comparing it with the 
standard RPA method in section IV. In particular, we 
demonstrate that we can use our scheme to study large 
size systems, beyond the capability of the standard RPA 
approach. Section V presents the results obtained by ap- 
plying the scheme to the model of Section II, including 
the density of states and nature of the eigenvectors of the 
magnetic excitations of the coupled fermion-spin system. 
We summarize our conclusions in section VI. 



II. THE MODEL 

In the following we develop the formalism for an im- 
purity band Hamiltonian that we believe to be appro- 
priate for describing (at least qualitatively) the diluted 
magnetic semiconductor Gai_xMnj;As for low doping x 
(in the insulating phase, or not too far from the metal- 
insulator transition). However, generalizations of this 
method for other types of Hamiltonians, crystal struc- 
tures, parameters etc. can be done in a straightforward 
manner. 

The III-V host semiconductor is assumed to have the 
zinc-blende structuie appropriate for GaAs. Experimen- 
tal measurementsEj suggest that valence-II Mn substi- 
tutes for valence-III Ga. As a result, the Nd Mn dopants 
are randomly distributed at positions Ri, i = 1 , . . . , 
on the N X N X N FGG Ga sublattice, corresponding 
to a; = Nd/4:N^. All throughout this paper, we assume 
periodic boundary conditions. 

Due to the valence mismatch, each Mn introduces a 
charge carrier (hole) in the system. However, experimen- 
tally it is found that the hole concentration is consider- 
ably suppressed through compensation processes. As a 
result, the number of holes is only a fraction p = 10 — 30% 
of the number of Mn atoms Nh — pNd- Each Mn atom 
also has a |— spin Si from its half-filled 3d level. The 
magnetic properties of the doped semiconductor are re- 
lated to the exchange interaction between the Mn spins 
and the charge carrier spins. This interaction is known 
to be antiferromagnetic (AFM), and proportional to the 
probability of finding the charge carrier at the corre- 
sponding Mn site. 

Recently, we proposed a simple impurity-band model 
to describe the behavior of the charge carriers for the 
low-doping regimeJ13't3 The main justification is that 
near and below the metal-insulator transition (x ~ 0.03) 
the density of charge carriers is not large enough to ef- 
fectively screen the attractive Coulomb potential of the 
Mn dopants. As a result, bound, hydrogen-like impu- 
rity states are created about each Mn site, at an energy 
Eh ~ 1 Ry above the top of the valence band. Due to in- 
teractions, these impurity states lead to the appearance 
of an impurity band, and the holes first occupy states 
in this band. Only if the concentration of holes (or the 
temperature) is large enough, are states in the valence 



band itself occupied by holes. Thus, it is reasonable to 
attempt a description of the charge carriers behavior only 
in terms of this impurity_hand. 

As in previous work,E£lE^ in the following we will use 
the electron formalism to treat the problem. This is 
equivalent to performing a particle-hole transformation 
which leads to an inversion E — > —E of the charge car- 
rier spectrum. Thus, instead of emptying the top of a 
valence-like impurity band (i.e., introducing holes in the 
system) we instead occupy the bottom of a conduction- 
like impurity band. We also make the simplifying as- 
sumption that the isolated impurity wave-function for a 
charge carrier trapped near the Mn at site Ri is the Is 
hydrogen orbital 0(r—i?i) = </)i(r) exp(— |r — Ri\/aB), 
where qb = h/ \/2m*Eb = 8 A is the appropriate 
Bohr radius related to the effective mass of the heavy 
hole (m* « 0.5me for GaAs) andr-the binding energy 
{Eb = 112.4 meV for Mn in GaAs).Ea Our approach ne- 
glects both the complicated orbital form of the acceptor 
wave-function (Ref. p6| ) and spin-orbit coupling. The for- 
mer is not expected to lead to any qualitative changes. It 
has recently been proposed that spin-orbii-coupling leads 
to frustration in the magnetic ordering.E3 These effects 
are left out in the present study, which concentrates on 
the effect of disorder. A study combining both effects is 
indicated for future work. 

We consider the Hamiltonian 

H = ^ tijcl^Cja- + ^ JijSi ■ Sj. (1) 

Here, creates a charge carrier with spin a in the impu- 
rity state centered at site Ri. The first term describes the 
hopping of charge carriers between impurity states. We 
use the parameterization tij = 2(1 -I- r/ae) exp (— r/ae) 
Ry, where r — \Ri — R^ l, of magnitude and form ap- 
propriate for hopping between two isolated Is impurities 
which are not too close to one another.cZl This hopping 
term has been shown to lead to the appearance of an 
impurity band which has a mobility edge, as well as a 
characteristic energy for thex)ccupied states in agreement 
with physical expect at ion .Es 

The second term of the Hamiltonian ([^) describes the 
AFM exchange between the Mn spin Si and the charge 
carrier spin Sj — \'^ja_'^a(iCjf3 [t? are the Pauli spin ma- 
trices] . This AFM exchange is proportional to the prob- 
ability of finding the charge carrier trapped at Rj near 
the Mn spin at Ri, and therefore Jy = J\(f)j{Ri)\^ — 
Jexp (^—2\Ri — Rj\/aB^. Based on calculations of the 
isolated Mn impurity in GaAs, we estimate the exchange 
coupling between aJaole-and the trapping Mn {Ri = Rj) 
to be J = 15 meVJl3'L3 As already stated, the number 
of Mn atoms is Nd, and therefore there are a total of 
2Nd states in the impurity band. The number of charge 
carriers is fixed to iV?, = pNd, where we take p — 10%. 

In Ref. |l^ we studied the relevance of various other 
terms, such as an on-site potential (associated with 
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the Coulomb potential of the compensation centers), 
Hubbard-like on-site repulsion (describing interactions 
between charge-carriers), external magnetic field etc. 
While they lead to various quantitative changes, we be- 
lieve that the qualitative picture we present in the follow- 
ing sections is not changed by their absence. 

III. THE COLLECTIVE MODE SPECTRUM 

The Random Phase Approximation (RPA) describes 
the collective excitations of a system about its self- 
consistent T = mean-field ground state. In order to 
clarify the notation used, we begin this section with a 
very short review of the derivation of the relevant equa- 
tions for the mean-field ground state. 

A. The mean-field ground state 

We use the equation-of-motion approach of Ref . |2^. At 
the mean-field level, we are trying to find non-interacting 
fermionic quasiparticles through a unitary transforma- 
tion of the on-site charge carrier operators: 

aL = E^»-Wcl (2) 

i 

such that 

[T-{-,aicr]^ Enaal^ (3) 

and 

[n,s+]^-hH,s+ (4) 

For a non-interacting system, Eqs. (||) and (^ would be 
exact. For an interacting system at the mean-field level, 
the exact Hamiltonian is approximated by the diagonal- 
ized, non-interacting Hartree-Fock Hamiltonian 

H HhF = ^ Enaai^ttna - ^ HiS- (5) 
na i 

for which Eqs. (||) and (^) become exact. 

It is straightforward to show that for Hamiltonian (|l|) 

[T-t-: aL] = Yl '^•?')'=^<^ ^j-J'^'"' 

ij ij 

The right-hand-side of Eq. (||) contains spin operators, so 
it cannot be put in the form of Eq. (^) unless we linearize 
it, by replacing the spin operators with their mean- field 
ground state expectation value Sj (Sj) ~ e.zSMn{j)- 
While the choice of a coUinear ground-state, with all Mn 
spins aligned along the 2;-axis, is not the most general 
possibility, we have actually shown in Ref. |l^ that the 
ground-state of this model, for the range of parameters 



we are interested in, is indeed coUinear. This is why we 
choose to apriori make this assumption in this case. 

After the linearization, we can now equate the right- 
hand-side of Eq. (H) with EnaO-na [sGC Eq. (^] to obtain 
the Hartree-Fock equation for the one-particle orbitals: 

Enatpna{i) = Ujlpnajj) + Jij S Mn{j)^1pnaii) (7) 

j j 

The mean-field ground-state of the charge carriers is ob- 
tained by occupying the Nh lowest-energy orbitals, 

l*>cc- n ''nalO) (8) 
{na} = l 

The effective magnetic fields Hi [see Eq. ||)] are ob- 
tained in a similar way. We evaluate the commutator 
of the exact Hamiltonian with the spin raising operator 

[W. = E J^^^^'t + J^^^t^t (9) 

The right-hand-side now contains charge carrier opera- 
tors (in the charge carrier spin operator) , beside the spin 
operators. We again use linearization sf — > {sf) = [see 
Eq. (^], si (sf) = s/i(i), and by comparison with Eq. 
d) we find 

H, = -Y,J^oSh{J) (10) 

3 

From Eq. (^ we see that if Hi > 0, then in the mean- 
field ground state the spin Si is in the "up" state | -I- S) 
and therefore SMn(i) — +'5', while if < the spin Si is 
in the "down" state | - S) and SMn{i) = -S [S ^ 5/2). 
As a result, the spin-part of the mean- field ground-state 
can also be easily found if the expectation value of the 
charge carrier spins Sh{j) are known from Eqs. (|^), (||). 
Thus, we obtain the usual self-consistent Hartree-Fock 
equations. 

Diluted Magnetic Semiconductors exhibit ferromag- 
netism at low temperatures. As a result, we assume that 
in the mean-field ground state, all Mn spins are fully po- 
larized SMn{i) = +S. Then, from Eq. (j^) we find the 
lowest-lying Nh states. If the charge carriers are also 
fully polarized at T = 0, all the occupied Nh state have 
(T =1, and therefore the mean- field ground-state is of the 
form: 

m = Y[ai^\Q)®\S,S,....,S) (11) 

n=l 

Of course, one must check that self-consistency is 
obeyed by verifying that indeed the first Nh lowest en- 
ergy one-electron states are all spin-down. We find that 
this always holds true for the parameters we study in 
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this paper. As discussed later, the spin-wave spectrum of 
the collective excitations also confirms that this ground- 
state is indeed stable. However, for higher charge carrier 
concentrations (or Hamiltonians with other types of in- 
teractions) it is likely that the HF ground-state is only 
partially spin polarized. In that case, one must do the 
full iterational search for the self-consistent HF ground 
state. The computation for the spin-wave spectrum in 
the partially-polarized case (for instance due to spin-orbit 
coupling), can be derived in a similar way to the one we 
present in the following for the fully polarized case. 



B. The Random Phase Approximation 

We will use the same equation-of-motion approach to 
derive the RPA equations. We are interested in spin- 
waves, which are related to spin-flip (spin-lowering) pro- 
cesses. Therefore, the RPA operators for spin- wave col- 
lective modes must be constructed in terms of spin-flip 
operators 




(12) 



The index h = 1, runs over the empty ("hole") states 
of the HF ground state, while the index p = 1, Nh runs 
over the occupied ("particle") states of the HF ground 
state. The index i = l,Nd runs over all the Mn spin 
positions. The HF ground-state j^*) is fully polarized, 
and given by Eq. (|l|). 

The spin-flip operators have "bosonic" nature, as ex- 
pected for collective modes operators. Indeed, it is 



straightforward to verify that 



— Spp'Shh' 



— Sij, with all the other commutators iden- 
tically zero. Here, the notation (...) signifies an average 
over the HF ground state (^'|...|^'). 

The Hamiltonian can be rewritten in terms of the a„o- 
operators as [see Eqs. (|l|), (0)] 



(13) 



where 



and 



tnma = ^ Co- («) V^m^r ( j) 



Jna,-m0{i) = ( j) (j) • (14) 



One can easily derive the equations of motion for the 
spin-fiip operators. For instance 



(O'^aLam.T'S', (15) 



Since RPA is an approximation, not an exact solution, 
we again must linearize this commutator about the HF 
ground-state, and keep only the meaningful terms. In 
other words, we replace 

Here, the first transformation is the linearization about 
the HF ground state | . The second transformation re- 
sults if one uses the HF expectation values {Sf) = S, 
(ajj|a„i|) = 0. The restriction to relevant terms means 
that in any sum over (n |), we restrict n to being an occu- 
pied orbital in the HF ground-state 1 < n < Nh, because 
otherwise = O7 i-S- a spin-flip process from the 

HF ground-state is impossible [see Eq. (1^)]. Performing 
a similar linearization for al^^ama-S^ and collecting the 
various terms, we find the linearized equation of motion 

[H, Bj] « H,Bl + V2SY,J2 JplMi^)bhp (16) 
p=i h=i 

After similar steps, the linearized equation of motion for 
the other spin-flip operator bhp is found to be 

[H, bHp] « - {Eh^ - Ep_i)bhp-V2SY, JtiM^^B] (17) 



Eqs. (|16|) and dlT]) allow us to "guess" the RPA part 
of the Hamiltonian, for which Eqs. (^6|) and ( p7| ) become 
exact, to be given by 

Urpa = Y.T. (^"T - Epi) blphnp + H.BjB, 

h—lp—l i 



iVrf Nh 

+ V2SYY.Y.^JpiM^Whp + h.c.) (18) 

i h=lp=l 

Thus, the RPA Hamiltonian is quartic in electron oper- 
ators and quadratic in spin operators, showing that it is 
the next order term after the Hartree-Fock Hamiltonian 
(which is quadratic in electron operators and linear in 
spin operators) 

Ti ~ 'Hhf + Ti-RPA + ■■■■ 

Higher order approximation terms describe interactions 
between the collective modes, which are neglected at the 
RPA level. 

We want to diagonalizc the RPA Hamiltonian to the 
canonical form 



'HrpA — ^ ^T-^aQaQa 



(19) 
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where Qj^ is the creation operator for a spin-wave with 
energy hflai and a is an index (in homogeneous sys- 
tems, it would be a wave- vector) . Since they create 
collective mode excitations, these operators must have 
bosonic nature [QonQp] — 0. However, since they do 
not describe exact, but only approximative solutions of 
the exact Hamiltonian, in fact only the weaker condition 
{[QcQIi]) — <5q/3 is obeyed. The most general form for 
these operators [see Eq. 



can also denote 



(20) 



h=i p=i 



Using the equation of motion [Hn^pA, Qn] = ^^aQa and 
computing the commutator using Eqs. (p^), (|2C|), wc find 
the diagonalization condition to be 

i 

Nh. 

- m r/") = -V2SJ2 E (22) 



= 1 h=l 



These are the RPA equations and they can be recast ia. 
the customary standard eigenvalue RPA equation formE3 



X 



E J 

J* H 



(23) 



where the vectors X and Y contain all the unknowns X^p 
and Yi, and the matrices E, H and J have the elements 

and 



{Eh] - Epi), Hi 



'^hpji'p' 

3hp,i = V^Jpi,h]{i)- 

The dimension of the RPA matrix (^3|), and there- 
fore the number of normal modes, is Nd + Nh x Nd- 
Of these solutions, Nd are proper spin-wave collective 
modes, while the rest of Nh x Nd are spin-flip processes 
associated with particle-hole excitations. If there were 
no interactions (J=0), the eigenenergies of these spin- 
flip processes would be Hi for the lowering of the Mn 
spin i, respectively Eh] — E„\ for a hole spin- flip [see the 
left-hand side of Eqs. (21), (E3)]. Interactions renormal- 



ize these values [as described by the right-hand-side of 
Eqs. (^l|),(22)], but one still expects Nd spin- wave col- 
lective solutions at low energies, coming from the Mn spin 
lowering, and Nh x Nd spin-flip particlc-hole excitations 
at energies comparable or larger than the spin-flip gap 
A = Ei^ - Em^i- 

We now comment on the sign of the RPA spectrum 
frequencies hria. Since Hi > [see Eq. ([lO|)] and 
Eh] — Epi > (by definition of the HE ground-state), 
it is apparent from Eqs. (|2l|),(p2|) that Nd spin- wave so- 
lutions will have positive energies hfla , while the Nh x Nd 
spin-flip particle-hole excitations will have negative ener- 
gies. This does not imply that the system is unstable, 
but simply that we have not chosen the proper spin-flip 
creation operators. Instead of the choice of Eq. (pO|), we 



-Qi-T.il ^>hp + T yt'Bi (24) 

h=l p=l i 

In terms of these new operators, we have now 

\HRPA,Pa\ = TT-^aPa, 1-6. 



Hrpa = ^(-?ir!a)PiPa 



(25) 



In other words, a negative solution for hfl simply means 
that we chose the corresponding annihilation operator 
instead of the creation operator when we wrote Eq. (|2^). 
For spin- flip processes, it is obvious that this is the case 
from the dependence of Qj^ on bhp, instead of bj^^. 

An unstable mean-field ground-state is sigBa.]^ by 
complex values of the spectrum frequencies Tifla E3'Ej Ex- 
cited states about the mean- field ground-state |^) are of 
the general form |$) = [l-f E Ca(Ql)""]|^')- Within 
RPA the dynamics of such states is dictated by Hrpa, 
leading to a time dependence: 



mt)) 



[1 + y^Cae" 



where Ehf = (^'|7Y|5') is the Hartree-Fock energy of the 
system. If any of the frequencies fla has a non-trivial 
imaginary part, it follows that the expectation value of 
any operator (<I>(t)|...|$(i)) will move in time exponen- 
tially away from its mean-field ground-state expectation 
value (\E'|...|5'), i.e. the mean- field is unstable to small 
perturbations. 

The advantage of the standard RPA approach is that 
by solving the eigenvalue problem ( p3| ) , one finds the nor- 
mal mode spectrum TiQ and the spatial distribution of 
the spin-waves (related to the Yi coefficients) at once. 
The obvious disadvantage is of a numerical nature: for a 
disordered system the RPA matrix must be diagonalized 
numerically. The size of the matrix is n = (1 -t- Nh)Nd. 
The typical sizes we consider are systems with around 
Nd = 500 Mn spins, and Nh = 50 = 10%Nd holes (al- 
though concentrations up to 30% might have to be con- 
sidered, depending on the doping x), leading to RPA 
matrices of typical size n > 25, 000. As a result, one has 
to either consider much smaller systems (in which case, 
finite size effects may be overwhelmingly important), or 
to try sparse matrix techniques to obtain the first few 
low-energy modes of the RPA matrix (although the off- 
diagonal J matrix is not necessarily sparse). 

There is, however, an alternative way of finding the 
collective mode spectrum and spatial distributions. We 
can directly solve for the Xhp coefficients from Eq. (|21 
and rewrite Eq. (|2^) as 



'^Mij{uj)Yj{Lu) = hLjY,{uj) 



(26) 
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where 



issues related to instability. 



p=i h=i 



TiLO + Eh] — Epi + irj 



(27) 



The advantage of this formulation is that one has to 
deal with much smaller matrices (the dimension of the 
M matrix is = 500). Eq. (|2^) tells us that for fre- 
quencies ri belonging to the spin- wave spectrum, the ma- 
trix M(ri) has at least one eigenvector corresponding to 
an eigenvalue \{^a) — fi^la (if the mode is degenerate, 
there are several such eigenvectors). The corresponding 
eigenvectors give us the desired spatial mode distribution 
y}"^ = Y.iria). Thus, the proble m is reduced to sweep- 
ing the range of frequencies of interest (for low-energy 
collective excitations, this is usually a fairly small range 
of frequencies near uj = 0), and monitoring the eigen- 
values of the M([x)) matrix. The dependence on uj of the 
matrix elements AIui {uj), and therefore of the eigenvalues 
X{lli) is monotonic for small a; <C A [see Eq. (p7|)]. As 
a result, each equation X{lu) — hut has a single solution, 
and we expect to find Nd collective mode eigenenergies, 
one for each eigenvalue. The monotonic dependence of 
X{uj) on u! also simplifies the search for the collective spec- 
trum, since it is enough to evaluate the eigenvalues X{uj) 
in the range of interest on a grid of step 6uj, and use 
linear interpolation to find the solutions of the equations 
X{lu) — huj. As we show in the following by comparison 
with standard RPA (for small system sizes) we can very 
easily achieve relative errors less than 10~^. 

The RPA spectrum has a second type of solutions, the 
spin-flip particle-hole continuum. These solutions are at 
frequencies of the order of the gap between the last occu- 
pied and first empty HF states TujJ ~ A, and one expects 
Nh X Nd such solutions. In principle, one can use the same 
method to find these energies as well, although the fact 
that they extend over a larger range of frequencies and 
that the uj dependence is very strong makes the search 
more difhcult. Typically, one is only interested in the 
lower and upper limits of the spin-flip continuum, which 
can be found easily with the method just described. 

Finally, we would like to mention that it is always pos- 
sible to reformulate the standard RPA equation (|2^) in a 
"self-consistent" form of the type shown in Eq. (26), 



with a matrix M(aj) of much smaller dimension than 
the RPA matrix. For instance, for Hamiltonians de- 
scribing interacting electron systems, the RPA matrix 
has a dimension ^ NgNe, where No is the number of 
occupied, and iVg is the nuHiber of empty orbitals in 
the mean-field ground-state.E^ For a system with a finite 
concentration of electrons x on a lattice of linear dimen- 
sion in a D-dimensional space, we have No ~ xN^ 
and iVe ~ (1 — x)N^, and the RPA matrix scales as 
x{l — x)N'^^ . For such systems, one can always find an 
equivalent reformulation of therB,PA equation in terms 
of a M((jj) matrix of size ^ N^rB This allows numerical 
handling of considerably larger systems, without having 
to resort to sparse matrix techniques, which often have 



IV. IMPLEMENTATION 

In this section we show, by direct comparison against 
known calculations, as well as by direct comparison 
against solving the RPA matrix equation, that the for- 
mulation of Eqs. (p6|), ( [2^ ) gives correct and numerically 
very accurate results. In the second part of the section 
we describe in more detail the numerical implementation 
as well as efhciency of our method. 



A. Analytical solution 

We first apply our approach by solving Eqs. (^, ( [27| ) 
analytically for a simplified case. We assume that the Mn 
spins are arranged on an ordered superlattice, instead of 
having random positions. Then, the charge-carrier part 
of the mean-field Hamiltonian is easily diagonalized in 
the fc-space: 



where 



ka 



Er = e,- H cr 

ka k 2 



(28) 



(29) 



Here, e^r = X^a^o ^5 ^xp (ik ■ 5) is the kinetic energy of the 
non-interacting electrons, where 5 indexes all the neigh- 
boring sites and t-g — tij for which \Ri — Rj\ = \d\. 
Also, Ja = "^(5' '^here again Jg ~ Jij for which 
\Ri — Rj\ = \5\. In the HF ground state, the states (fc, |) 
given by the condition Nh = J2\k\<kF ^ occupied. 
Since the Nh electrons are equally distributed among the 
Nd Mn sites, and are fully polarized, the z-component of 
the total spin created by electrons at each site, within 
HFA, is Sh{j) -- , _ 
H, = Jop/2 [see Eq. (0)] 

Using the fact that the solutions of Eq. (^6| ) must 
also have plane- wave structure Yi ~ exp (iQ ■ Ri), it fol- 
lows that the index a Q becomes a wave- vector, as 
expected, and we can reduce Eqs. (|2^) and (|2^) to a 
self-consistent equation for each collective mode fiLo{Q): 



= — p/2, and therefore we find 



flUJ^ 



2 



S\JiQ)\' 

2Nd 



k 



^k-Q 



JoS 



(30) 



Here, J{Q) = ^^exp (iQ • (5) J^. The occupation num- 
ber fiEj: ^) obviously comes from the sum over occupied 
states we had in Eq. (p7|). For a finite-size ordered lat- 
tice, one expects considerable degeneracy for each mode 
Ej:^, due to invariance to various symmetry transforma- 
tions of the fc-vector. In most cases, the number of charge 



7 



carriers is such that the last orbital (of degeneracy G) is 
only partially occupied by g charge carriers. In this case, 
one must obviously choose / = g/G for all these orbitals, 
and f = I for lower, fully occupied orbitals, and f — 
for higher, empty orbitals. Otherwise, the translational 
invariance is broken and the collective spectrum is not 
indexed by a wave- vector. 

The sum in Eq. (|3^) can be performed if we as- 
sume that the dispersion at the bottom of the band is 
of quadratic form eg = h^k'^/{2m) (this is a reasonable 
approximation since we are interested in low filling frac- 
tion p). Then, the sum over occupied states can be trans- 
formed to an integral which is straightforward to evalu- 
ate, leading to the solution: 



^VfQ 



(31) 



Here, vp = hkp/m is the Fermi velocity, where the Fermi 
vector is given by 



Nh - pNd 



(27r)3 3 ■ 
The function f{Lo^ Q) is given by 

hui + €Q + JqS 



huj + eg + JqS 
vpQ 



VfQ 



In 



huj + eq + JqS — vpQ 

hu; + €q + JqS + VpQ 



where eg = h^Q^/2m. 

For comparison purposes, we will assume a local AFM 
interaction cJpdSij, leading to J{Q) — Jq = cJpd- 

In this case, the Hamiltonian becomes identical to the 
Hamiltonian used by Konig. et. al in Ref. ^ provided 
that the effective mass m in the dispersion relation is 
assumed to be the band effective mass. We can easily 
solve Eq. (|T]) in the asymptotic limit Q ^ 0, to find 



hioiQ) 



1 



2m 2S/p - 1 



P 



4 t 
^~5J^S V67r2 



2/3' 



(32) 



This is indeed identical with the asymptotic long- 
wavelength spin-wave spectrum obtained in Ref. ^ and 
has the typical quadratic dependence of spin-waves in 
conventional ferromagnetic systems. 



B. Direct comparison between the two 
formulations 

In this section we briefly illustrate the accuracy and 
speed of our formulation of the RPA problem [Eqs. ( |2^ ) 
and (2^)], in comparison with the standard RPA ap- 
proach[Eq. (p3|)]. We use a rather small system, with 
Nd = 80 Mn spins randomly distributed on a FCC lattice 




6.0 

to (meV) 



FIG. 1: The largest eight eigenvalues X{uj) of the M(lj) 
matrix (circles), evaluated on a grid with a step Suj = 0.05 
meV. The full line is fiiD. The spin- wave spectrum is given by 
the condition \uj — Ulo. 



with linear size A'^ = 10, corresponding to a Mn concen- 
tration X = Nd/4:N^ = 0.02. The number of charge 
carriers is — pNd = 8, and the other parameters are 
as defined in Section II. Thus, the standard RPA in- 
volves the diagonalization of a non-hcrmitian matrix of 
dimension 720. 

In Fig. we show the w-dependence of the largest 
eight eigenvalues A(a;) of the M matrix (circles), evalu- 
ated on a grid with a step Suj — 0.05 meV. [For real oj and 
r] = 0, the matrix M(w)is hermitian and all eigenvalues 
A(ti;) are real, see Eq. (p7|)]. The full line is huj, and the 
spin-wave spectrum is given by the condition X{uj) — huj. 
One can see that the eigenvalues A(ix') have monotoni- 
cally increasing dependence on oj for small w ^ A « 50 
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FIG. 2: Accuracy Alu/uirpa for all Nd = 80 spin- wave fre- 
quencies ujhpa, between the exact RPA (standard) values and 
values obtained in our formulation using a grid with Sco — 0.05 
[see Fig. 
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FIG. 3: Largest relative error /S.lu/ujrpa among 5Nd ~ 400 
spin- wave frequencies ujrpa, between the exact RPA (stan- 
dard) values and values obtained in our formulation using a 
grid with 5lo = 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.75 and 1 meV. The 
line is drawn to guide the eye. Even for a very large grid, the 
relative error is still reasonably small. 



meV, and therefore the equation X{uj) — Tilo can have 
at most one solution for each eigenvalue, for small u. 
If each eigenvalue yields a solution, we have found all 
Nd spin- wave frequencies. If one or more eigenvalues do 
not intersect the fiuj line, this means that the mean-field 
ground-state is unstable ( there are complex spin-wave 
frequencies). For all the cases and parameters we inves- 
tigate, we found that the fully-polarized ground-state is 
stable for our model. 

To find the spectrum frequencies ?if2Q, we use linear 
interpolation over the interval 5lo where each A(w) curve 
intersects the Tilo curve. This method avoids the need 
of evaluating the eigenvalues at too many uj points (the 
most time-consuming part of the computation is the eval- 
uation of the matrix elements of M) . The high degree of 
accuracy obtained with this method is demonstrated in 
Fig. ^, where we compare the exact RPA spectrum ujrpa 
obtained through direct diagonalization of the RPA ma- 
trix [see Eq. (Pq)] with the the values Cjrpa obtained 
through our formulation. In fact, we plot the relative 
error Au/ujrpa = (Corpa ~ ujrpa)/i^rpa for each spec- 
trum frequency ujrpa- As one can see from Fig. ||, the 
largest error is of the order 2-10^^, for the grid Slo — 0.05. 

This suggests that one could use a much larger grid 
5lj and still have a good accuracy. Indeed, we have com- 
puted the spectrum through both methods for five dif- 
ferent realization of disorder, using various grid values 
5liJ in our formulation, and selected the largest relative 
error for each case (from a total of bNd — 400 values). 
The results are shown in Fig. ^. The error is found to 
scale as the square of the grid size. Even with a very 
large grid Slo — \ meV, (which implies the evaluation of 
the eigenvalues A(w) in only very few points), we still get 
a relative error smaller than 10"'^. Thus, one can use 
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FIG. 4: CPU time for the two RPA methods, for systems 
with Nd = 64, 125 and 216 spins. The standard RPA (full 
line) is significantly more time consuming than our formula- 
tion (dashed line). 



the grid step 5lj to optimize and considerably speed up 
the computation. However, as the number of spins (and 
eigenvalues) Nd increases, one must take into considera- 
tion other complications, as discussed below. 

In Fig. ^ we compare CPU times for finding the RPA 
spin-wave spectrum using the standard vs. our RPA for- 
mulations. All simulations were done on the same pro- 
cessor. We used a small grid 5bj = 0.05 meV for our 
approach, and verified that all relative errors were less 
than 10"^. We use systems with Nd = 64, 125 and 216 
spins, randomly distributed on FCC lattices of linear size 
N = 12, 15 and 18 [x = 0.0092, p = 10%). As expected, 
the computational time depends in a power-law manner 
on the size Nd, with exponents 6.4 and 3.2 respectively 
for the standard and current approach. Clearly, our for- 
mulation can be successfully used for much larger sizes 
than those afforded by the standard approach. 

One important aspect to keep in mind is that as the 
system size Nd increases, neighboring eigenvalues become 
more closely spaced, and accidental crossings and anti- 
crossings appear. Some typical examples are shown in 
Fig. 1^, for a disordered system with Nd = 512 spins. 
Tracking the crossing of the eigenvalues is essential, since 
such events lead to a change in the indexing of the eigen- 
values from one grid point to the next. We found that re- 
quiring continuity in first and second derivatives allows a 
unique identification of each continuous eigenvalue from 
one grid point to the next (3 steps are shown for each 
eigenvalue in Fig. ^). In fact, the only relevant crossings 
are the ones that take place within the step 5uj where 
the eigenvalues intersects the huj line (shown as a dotted 
line in Fig. m. Such an example is shown in the left 
panel of Fig. ^. For the largest system size investigated, 
of Nd — 512 randomly distributed spins, we find that 
only 4.5% of the eigenvalues have such relevant cross- 
ings within Slo of A(u;) = hu. The percentage decreases 



9 



0.18 




0.14 



0.138'-; 



0.136- 



0.164, 



0.2 0.4 0.6 O.i 

0) (meV) 




1 0.2 0.4 0.6 o.i 

at (meV) 



0.8 



„0.6 

> 

<1> 
E, 

•5^0.4 

a 

0.2 



o 


N = 


15 


□ 


N = 


18 


A 


N = 


24 




N = 


150 




(-71,0,0) (11,0,0) (-11,-71,0) (Jt,7t,0) (-71,-71,-71) 

(j-space 




(71,71,71) 



FIG. 5: Eigenvalues A(tj) for a system with Nd ~ 512 random 
spins, evaluated for a grid SLu—0.05meV . Crossing and anti- 
crossings of neighboring eigenvalues is apparent. The dotted 
line is hu. Since we are interested in the solutions of the 
equations X{uj) = hui, it is clear that too large a grid 5uj may 
lead to considerable errors . 



with decreasing Nd, to 1.1% and 0.3% for Nd — 216 and 
125, respectively. This percentage also depends on the 
grid step Slu; for increasing Slo the identification of cross- 
ings and anti-crossings becomes more difficult, leading to 
possibly large errors. One can optimize the choice of the 
grid step Slo by starting with a larger value. The well 
separated eigenvalues (such as the ones depicted in Fig. 
|l|) will provide unique identification and very accurate 
values for their corresponding spin- wave spectra values. 
However, where considerable mixing and therefore non- 
linear variation of the eigenvalues A(a;) is apparent, a 
finer mesh is necessary in order to correctly characterize 
their variation near hoj. In all our simulations, we use 
the grid Slo = 0.05 meV, which allows for comfortable 
tracking of each eigenvalue and also is sufficiently small 
to allow us to approximate the variation of A(a;) as being 
linear within each Slo step. 

These comparisons clearly demonstrate the accuracy 
and speed of our formulation of RPA, as compared to 
the standard RPA. The biggest advantage, though, is 
that it can easily be applied to systems with large sizes, 
for which standard RPA is numerically cumbersome. 



V. SPIN- WAVE SPECTRA OF DMS 

In this section, we use the RPA method described 
above to study the density of states (DOS), as well as 
nature (extended or localized) of the spin- wave spectrum 
of (III,Mn)V DMS described by the model introduced in 
Section II. 



FIG. 6: Spin-wave dispersion Q{q) for an ordered, simple 
cubic superlattice arrangement of Mn spins inside the host 
semiconductor. Lattices of linear size A'^ = 15, 18, 24 and 
150, with a total of respectively 125, 216, 512 and 125000 Mn 
spins are considered. This corresponds to x = 0.0092 and 
p = 10%. Dispersion is plotted along three linear cuts in the 
Brillouin zone: (— 7r,0,0) (7r,0,0), (— vr, — 7r,0) (7r,7r,0) 
and (— TT, — TT, — 7r) (7r,7r,7r). 



A. Spin-wave density of states 

We study the spin-wave density of states for three types 
of arrangements of the Mn spins inside the host semicon- 
ductor. All samples studied correspond to a; = 0.0092 
and p = 10%. Other parameters are as specified in Sec- 
tion 11. 

First, we consider fully ordered systems, in which the 
Mn spins are arranged on a simple cubic superlattice. For 
X = 0.0092, the superlattice constant is ql = 3a. In this 
case, we can study the spin-wave dispersion and density 
of states for very large sizes, since we can use directly 
Eq. (|3^) to compute the spin- wave frequency fiuj{q) for 
each wave- vector q inside the Brillouin zone. Spin- wave 
dispersion obtained along three cuts in the Brillouin zone 
is shown in Fig. ^. The dispersion has quadratic behavior 
near the center of the Brillouin zone, as expected from 
the discussion for an ordered case provided above. The 
finite-size effects are reasonably small. For the small sizes 
we used both Eq. (|3l| ) and our method to compute the 
dispersion. The results of the two agree with a relative 
error of less than 10^^. 

One important aspect to notice is the small range of the 
spin-wave spectrum, as compared to the AFM exchange 
J = 15 meV. This is a consequence of the fact that the 
Mn spins do not interact directly with one another. In- 
stead, their interaction is mediated by the rather small 
concentration of charge carriers present. 

We compute the density of states (DOS) p{E) asso- 
ciated with the spin-wave dispersion for the superlat- 
tice case employing the standard method of dividing the 
Brillouin zone into tetrahedra and linearizing the disper- 
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FIG. 7: Upper panel: average density of states p(logj^Q E) 
on a logarithmic scale for systems with Nd = 125, 216 and 
512 Mn spins in moderately disordered configurations (full 
lines). Lower panel: average density of states for systems 
with Nd = 125, 216 and 512 Mn spins in strongly disordered 
configurations (full lines). The dotted line is the spin- wave 
density of states of a DMS with fully ordered (superlattice) 
configuration of Mn spins. All samples correspond to x = 
0.00924, p = 10%. 



sion (Ref. 30). The DOS obtained for the lattice with 
N = 150, Nd = 125,000 is shown as a dotted line in 
Fig. We use a logarithmic scale for the energy, and 
normalize the density of states such that 



dxp{x) 



where x — \ogiQ{E). This convention will be maintained 
throughout the rest of the paper. 

In the upper panel of Fig. ^ we also plot the den- 
sities of states obtained for moderately disordered con- 
figurations with Nd = 125, 215 and 512 Mn spins (full 
lines). These are configurations in which we place the 
Mn spins randomly on the FCC Ga sublattice of the 
host semiconductor, subject to the constraint that the 
distance between any two spins is larger than 2a. [In the 
ordered cubic superlattice, nearest neighbor spin sepa- 
ration is ol = 3a]. This moderate amount of disorder 
breaks translational invariance, and the wave- vectors are 
no longer good quantum numbers. Also, the large de- 
generacies of the superlattice spectrum are lifted. We 



computed the spin- wave spectrum for 200 realizations of 
the disorder with Nd = 125 spins, 100 realizations with 
Nd = 216 spins and 50 realizations with Nd = 512 spins. 
Thus, we have a total of over 21,000 spin-wave energies 
for each size, and statistics can be comfortably carried 
out. In particular, from Fig. ^ we see that the DOS 
histograms corresponding to the three sizes are smooth 
functions, i.e. the average over disorder is well accounted 
for. Also, the curves are practically indistinguishable 
from one another, implying that finite size effects are 
negHgible. 

In the lower panel of Fig. |^ we plot the DOS corre- 
sponding to fully disordered configurations with Nd = 
125,215 and 512 Mn spins (full lines). These are con- 
figurations in which we place the Mn spins randomly on 
the FCC Ga sublattice of the host semiconductor, with 
no restrictions (except that each Mn occupies a different 
site) . The number of samples analyzed is the same as in 
the previous case. Again, the curves show smooth be- 
havior and no finite-size effects. The origin of the peak 
appearing near w ~ 5 meV will be discussed later. 

Before continuing the analysis, we must also point out 
that the Goldstone modes have been left out in the DOS 
shown in Fig. ^ The RPA system always has one solu- 
tion of energy lu — (numerically, we find its magnitude 
to be less that lO^^''), corresponding to Yi = 1/^/Nd- 
This is the Goldstone mode, describing the same over- 
all rotation of all the spins. These Goldstone modes 
would appear in the DOS as a 6{uj) function at lu = 
(logio^; = -oo). 

The effect of disorder in the positions of the Mn is re- 
fiected in the considerable widening of the DOS (on a 
logarithmic scale), and rounding of the van Hove singu- 
larities of the superlattice DOS, as the amount of disorder 
increases. More significantly, however, is the substantial 
enhancement of the DOS at low energies. This behavior 
is in agreement with the general expectation of the effects 
of disorder on an energy-band dispersion. 

A question of considerable interest concerns the nature 
of the spin-wave excitations. We know that the ordered 
superlattice has extended, plane-wave type spin-waves. 
The general expectation is that disorder will lead to lo- 
calization, and thus one expects the appearance of local- 
ized spin-waves as the disorder increases. One way to 
characterize the nature of the spin-waves is to compute 
their Inverse Participation Ratio (IPR), defined as: 



IPR„ = 



(33) 



For an extended mode a (plane- wave, for instance), we 
expect that all f/"^ coefficients are roughly of equal mag- 
nitude, since all spins are expected to participate equally 
in the spin-wave. Then, it follows that for an extended 
mode a, IPRq ~ ^/Nd, i.e. it is inversely proportional to 
the size of the system. On the other hand, in a localized 
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FIG. 8: Upper panel: average IPR(E) for systems with 
Nd = 125,216 and 512 Mn spins in moderately disordered 
configurations (circles, square and triangles, respectively). 
Lower panel: same for strongly disordered configurations. All 
samples correspond to x = 0.00924, p — 10%. 



mode a, only the spins within the locahzation volume 
have non- vanishing values for y/"^. Thus, it follows that 
for a localized mode a, IPRq is independent of Nd- In- 
stead, its value is given by the inverse number of Mn 
spins participating in the localized mode. 

We compute the IPR of the spin-waves in the follow- 
ing manner. As we sweep the frequencies lu of interest 
and compute the eigenvalues A(a;) of the M(w) matrix on 
the grid Stu, we also compute the corresponding IPR of 
the eigenvectors Y(w) (which are normalized to unity). 
The IPR is a continuous function of w, and we use linear 
interpolation to find its value at the spin-wave frequen- 
cies ritj of interest. Whenever eigenvalues X{uj) cross, 
the IPRs are no longer well defined: one can choose any 
linear combination of the eigenvectors of the degenerate 
modes, which would lead to different values for the cor- 
responding IPRs. We have checked that these acciden- 
tal crossings are rather rare (below a few percent) and 
equally distributed over the entire energy scale. If we ex- 
clude all these degenerate cases, we get a DOS which is 
virtually indistinguishable from the one obtained when 
these modes are kept. More importantly, we find that 
modes which cross predominantly have the same nature 
(either localized or extended). As a result, the ensemble 
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FIG. 9: Histogram of the IPR values for all spin- waves (spin 
flips) of energy E <5- 10"^), for Nd = 125, 216 and 512. Re- 
sults are shown only for strong disorder configurations. Ex- 
cept for the low-lying Goldstone modes, the histograms are 
identical, suggesting only localized modes at these energies. 



averaged IPR is not sensitive to the precise treatment of 
these cases. 

In Fig. H we plot the (geometric) average IPR(E) for 
the moderate (upper panel) and strong (lower panel) dis- 
order configurations. Again, the Goldstone modes are 
not shown. We use the geometric mean (i.e., arithmetic 
mean of the log values) in order to insure proper weight 
for the extended modes, with low IPR. For moderate dis- 
order configurations, we see that the spin-wave modes at 
high energies E > 1 meV are localized: the values cor- 
responding to the three different sizes collapse on top of 
one another. This is expected, given the fact that the 
upper band-edge of the superlattice spectra is just below 
E — I meV (see Fig. Thus, states of higher energy 
have been split off the band by disorder, and are ex- 
pected to be localized. The low-energy part of the band 
also appears to be localized: while the three curves do 
not quite meet, the value of the IPR is just below unity, 
showing that these spin-waves have considerable weight 
on only one spin (i.e., they correspond to a individual 
spin-flip). However, for the central part of the spectrum, 
the spin-waves are extended, with the IPRs for the three 
sizes clearly distinct and decreasing with increasing Nd- 

For the strong disorder configurations, this tendency is 
even more apparent. The high- and low-energy regions, 
E > 1 meV, respectively E < 5- 10~^ meV, contain local- 
ized spin-waves: in both limits, the IPR curves collapse 
on top of one another. The low-energy localized modes 
arc individual spin-flips. The associated Yi are vanish- 
ingly small at all sites except one, leading to IPR « 1. 
These sites are always situated far apart from all other 
Mn spins, and the probability to be visited by charge 
carriers is exponentially small (tailing from far occupied 
regions). As a result, the Mn spins at these sites are 
virtually isolated, and their spin excitations are Individ- 
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FIG. 10: Histogram of the IPR values for all spin- waves (spin 
flips) of energy E > 10"^), for Nd = 125, 216 and 512. Results 
are shown only for strong disorder configurations. The his- 
tograms are almost identical, suggesting only localized modes 
at these energies. 



ual spin-flips. The energy for such a spin-flip is equal to 
Hi , if one neglects small corrections due to the extremely 
weak interactions [see Eq. (p^ ) and following discussion] . 

Histograms of the IPR of all the modes with energies 
below E < 5 ■ 10^^) are shown in Fig. ||. Here, we 
also show the Goldstone modes, which have zero energy, 
and IPR = l/Nj,. The histograms have been scaled by 
the total number of modes for all disorder realizations 
considered for each particular size. Since the number 
of Goldstone modes exactly equals the number of differ- 
ent realizations of disorder considered, their peaks are in 
a ratio of roughly 200/100/50 with respect to one an- 
other. The main observation is that the histograms for 
the three sizes are very similar, with a huge peak just be- 
low IPR=1. The absence of a dependence on Nd confirms 
that all these modes are localized (except the Goldstone 
modes, of course). 

The high-energy localized states are of very different 
nature: their IPR is close to 0.5, suggesting that it has 
large contributions from only two sites. Analysis of the 
Yi values for such modes shows that most of them are 
associated with nearest-neighbor (n.n.) Mn spins on the 
FCC Ga sublattice. The characteristic energy for the 
spin-waves centered on such n.n. sites is just below 5 
meV, and this is the feature responsible for the peak ap- 
pearing in the DOS at these energies (see Fig. 0, lower 
panel). Careful inspection of the DOS reveals the ap- 
pearance of smaller peaks at somewhat lower energies 
(even the DOS for moderate disorder has a small peak 
at around 2 meV, and its IPR for these modes is close 
to 0.5). These peaks are associated with excitations of 
spin-pairs (or larger clusters) with varying separations, 
but they are not as well defined as the one corresponding 
to nearest-neighbors spins. 

The small clusters of Mn spins giving rise to such 
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FIG. 11: Histogram of the IPR values for all spin- waves (spin 
flips) of energy 5 ■ 10"^) <E < 10"^), for Na = 125 (upper), 
216 (middle) and 512 (lower panel). Results are shown only 
for strong disorder configurations. With increasing Nd the 
histograms shift to lower values, suggesting that these spin- 
wave modes are extended. 



modes are always found to be in regions densely pop- 
ulated with charge carriers. Due to their closeness in 
space, these spins are much more strongly coupled to 
one another than they are to other spins with which 
they share charge carriers. This leads to the "resonance- 
like" character of these modes: either of the spins can be 
fiipped with equal probability. As a result of the strong- 
coupling, the cost of flipping either spin is high, since 
it frustrates their FM arrangement. Indeed, the char- 
acteristic energy of roughly J/3 reflects much stronger 
coupling than the average one present in the superlattice 
case. The histogram of the IPR of all modes of energy 
E > 10~^) shown in Fig. |l^ conflrms all these conclu- 
sions. We caution that these energies may be substan- 
tially modified due to direct antiferromagnetic Mn-Mn 
exchange (left out in our model) in real systems. 

Finally, the spin-waves at intermediate energies 5 • 
10"'^ < E < 10~^ correspond to extended modes. It is 
rather difficult to establish exactly the "mobility edge" 
corresponding to the transitions from localized to ex- 
tended modes at either end, given the rather small system 
sizes considered here. Thus, we use the values quoted 
above only as plausible estimates of these boundaries. 
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Clearly, for these energies the average IPR decreases with 
increasing iV^i (see Fig. ||). In Fig. |l| we plot the 
histograms of the IPR for all the modes with energy 
5 • 10"'^ < E < 10~^. In this case, size dependence is 
clearly visible, both in the position of the broad peak 
maximum as well as in the lower cut-off in IPR value. 
The position of the maxima scale roughly like 1/Nd- We 
have attempted both gaussian and lorentzian fits, but 
neither seems to capture the low-IPR tail properly. In 
the smaller samples, a second peak appears just below 
IPR=1. At first sight, this might seem to indicate the 
existence of a finite density of localized states at these 
energies as well. However, we believe that this peak is 
a finite-size effect. Clearly, its amplitude decreases with 
increasing Nd and would vanish in the thermodynamic 
limit. This is very different from the behavior observed 
in either of the two localized energy ranges, where the 
distributions for all sizes are virtually identical (see Figs. 
^ and . The reason for the appearance of this second 
peak is simply the restriction IPR < 1. As the broad 
peak moves towards higher values with decreasing N^, 
all the values in the upper tail "bunch" at IPR=1. 

One final important observation relates to the abso- 
lute values of the IPRs in the extended spin-wave regime. 
Although scaling with system size is present, the corre- 
sponding IPR values are much higher than the ones ex- 
pected for modes extended over all Mn sites (these val- 
ues are shown by the positions of the Goldstone modes 
in Fig. ^, and are well below the cut-offs observed in 
Fig. ^l|). This means that in the disordered system, the 
extended spin-waves are delocalized only over a fraction 
of the total number of sites. As the amount of disorder 
increases, this fraction decreases, since IPR averages for 
the moderate disorder are smaller than IPR averages for 
full disorder. 

These results are consistent with the picture provided 
by the temperature dependent mean-field study of this 
model (Ref. 13). There, we concluded that in the dis- 
ordered system, the (small fraction p) of charge carriers 
are concentrated in a small volume of the sample, where 
the local Mn concentration is larger than average. As a 
result, the concentration of charge carriers is strongly en- 
hanced in these regions, and the exchange with the Mn 
spins inside these regions (which is proportional to the 
probability of finding charge carriers nearby) is greatly 
increased, preserving magnetization of these regions up 
to high temperatures. On the other hand, the Mn spins 
in the regions devoid of charge carriers are very weakly 
coupled, and behave as free spins down to very low tem- 
peratures. Due to the very inhomogeneous assembly of 
weakly- and strongly- interacting spins, the magnetiza- 
tion curves have very unusual, concave shapes. 

The present study of the spin-waves corroborates the 
same picture. We find very low-energy, spin-flip like ex- 
citations, which are obviously due to the weakly inter- 
acting spins, and which are responsible for a very sharp 
decrease of the, magnetization at exponentially small 
temperatures.E3ll3 This is is to be contrasted with the 



behavior of ordered conventional ferromagnets, where at 
low temperatures only low-energy long-wavelength spin- 
waves can be excited. Since their phase space is van- 
ishing in the long-wavelength limit (^ k'^ in d dimen- 
sions), the magnetization of conventional ferromagnets 
decreases very slowly from its saturation value with in- 
creasing temperature, leading to convex upwards magne- 
tization curves .EiJ 

In our model, the extended spin-waves are concen- 
trated around the high-density Mn regions, where the 
charge carriers mediating the interactions are to be 
found. This is consistent with the appearance of modes 
whose IPR, while scaling with the system size, shows 
modes extended only over a fraction of all system sites. 
As the amount of disorder decreases all the way to a fully 
ordered superlattice, the charge carriers are more and 
more homogeneously spread throughout the sample and 
the IPR of the extended modes has lower and lower val- 
ues, as observed from Fig. ^ However, the more homo- 
geneous a sample is, the less the average AFM coupling 
between the Mn spins and the charge carriers, since the 
charge carriers have now equal probability of being found 
anywhere in the sample, instead of being concentrated in 
a small fraction of the space. This leads to a decrease of 
the pritical temperaturft-Tc, as observed in both mean- 
fieldlla and Monte-Carldla analysis. The enhancement of 
the mean-field Tc with increasing disorder is suggested 
also by the existence of the high-energy localized modes, 
which show that ferromagnetic alignment will persist in 
high-density clusters up to very high temperatures. It 
also suggests that local ferromagnetic fiuctuations might 
be observed well above Tc- 



VI. FINAL REMARKS 

The aim of this paper is two-fold. First, we demon- 
strate the accuracy and speed of a new scheme for com- 
puting RPA spectra, which allows tackling of systems of 
much larger sizes than the ones that can be analyzed with 
the standard RPA formulation. Investigation of large sys- 
tem sizes and averages over many disorder realizations 
facilitate a clear picture for the problem we are inter- 
ested in, namely the spectrum and nature of spin-waves 
of a disordered DMS. 

We then demonstrate that disorder can significantly 
change the spectrum and the nature of the spin-waves. 
This is likely to lead to important consequences not only 
as far as magnetic properties are concerned (we have al- 
ready commented on the fast de-magnetization with in- 
creasing temperature, due to low-energy spin- flip modes) . 
More importantly, this may have significant consequences 
for transport properties as well. For instance, charge car- 
rier spin scattering is likely to be very different in vari- 
ous regions of the sample. LargCr-anomalous Hall effect 
has been observed in (Ga,Mn)As,CJ but the theory used 
to interpret it is borrowed from phenomenology relevant 
to homogeneous ferromagnetic metals. In an inhomoge- 
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neous system, some of the accepted ideas might have to 
be modified or at least verified to still hold true. 
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